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Abstract 

We evaluate the corrections to the matching coefficient of the vector current 
between Quantum Chromodynamics (QCD) and Non-Relativistic QCD (NRQCD) 
to three- loop order containing a closed heavy-fermion loop. The result constitutes 
a building block both for the bottom- and top-quark system at threshold. Strong 
emphasis is put on our completely automated approach of the calculation including 
the generation of the Feynman diagrams, the identification of the topologies, the 
reduction to master integrals and the automated numerical computation of the 
latter. 

PACS numbers: 12.38.Bx 12.38.-t 14.65.Ha 



1 Introduction 



A major goal of a future international linear collider (ILC) is the precise measurement of 
the top-quark production cross section close to threshold. Next to a precise extraction of 
the strong coupling, an unrivalled determination of the top-quark mass and its width is 
possible. This would open up a new chapter in the electroweak precision physics which 
leads to very strong checks of the Standard Model or possible extensions. 

The theoretical calculation of the threshold cross section is based on an effective the- 
ory [1,2] (for a review see [3]) which is constructed from QCD by integrating out the hard 
scale given by the top-quark mass. The connection between the two theories is established 
by so-called matching coefficients which constitute the coupling constants of the effective 
operators within NRQCD. 

A preliminary analysis to next-to-next-to-next-to leading order (NNNLO) of the top- 
quark threshold cross section, which is necessary in order to match the expected ex- 
perimental precision [4], has been performed in Ref. [5]. However, the three-loop static 
potential (see Ref. [6] for the fermionic contribution) and the three-loop matching co- 
efficient beyond the light-fermion approximation [7] are still missing. In this paper we 
provide a further building block needed for the completion of the NNNLO calculation: 
the heavy- fermion contribution to the three- loop matching coefficient. Next to important 
applications in the top-quark sector NRQCD is also an appropriate tool for the description 
of boundstate phenomena of charm and bottom quarks [3]. 



The matching coefficient of the vector current in the full and effective theory is defined 
through (k = 1, 2, 3) 



Q denotes a generic heavy quark with mass tuq and and x are two-component Pauli 
spinors for quark and anti-quark, respectively. In this paper we compute the three-loop 
non-singlet contribution to c v which contains one or two closed heavy quark loops. The 
analog corrections involving closed light (massless) quark loops have been considered in 



The evaluation of the three-loop diagrams contributing to c v is quite involved. One 
has to consider vertex diagrams with massive quarks on their mass shell and the external 
momentum q 2 = Atuq. This kinematical configuration in combination with an involved 
reduction to master integrals makes the calculation quite challenging. In this paper we 
discuss an automated setup which minimizes the manual interaction. Even the results for 
the master integrals are obtained in an automated way. 

Our approach for the automated calculation is introduced in the next section. Af- 
terwards we present the results for the matching coefficient in Section [3] and conclude in 
Section [U 




(1) 




(2) 



Ref. [7]. 
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2 Automated multi-loop calculation 



When evaluating multi-loop Feynman integrals one encounters several difficulties which 
have to be overcome. Among them are the generation of the Feynman diagrams, the 
reduction of the many integrals which appear at the initial stage of the calculation to a 
relatively small set of so-called master integrals and the evaluation of the latter. Very 
often the individual steps are automated, however, the interplay between them is not. In 
the following we present a setup where various program packages are combined in order 
to minimize the manual work. 

The individual steps of our automated setup are as follows 

1. All Feynman diagrams are generated with QGRAF [8] which requires two input files: 
one specifying the process and one containing the propagators and vertices occurring 
in the theory. 

2. The output of QGRAF is transformed to FORM [9] with the help of the program q2e [10, 
11]. q2e requires as input the FORM notation for the propagators and vertices and 
information about the hierarchy of the particle masses and the external momenta. 

3. The output of q2e is further processed with exp which identifies the various topolo- 
gies and generates for each diagram a separate file containing complete information 
like the projectors to be applied, the expansion to be performed and the topology 
file to be called. 

4. In a next step the FORM part is initiated. After taking the traces and applying the 
projectors a topology-specific file is included which expresses the result in terms of 
a sum of scalar integrals. 

5. From the sum of all diagrams a list of integrals is extracted. This list serves as input 
for crusher [12] which produces for each topology a table containing the reduction 
to master integrals. 

6. Each topology produces a certain number of master integrals. A small Mathematica 
routine combines all master integrals, identifies identical ones and generates relations 
among them. 

7. The tables from step 5 together with the relations among master integrals (step 6) 
are applied to the sum of the bare diagrams leading to a representation of the result 
in terms of a minimal set of master integrals multiplied with e-dependent coefficients. 

8. The input for crusher can also be used for FIESTA [13] which we employ in order 
to obtain numerical results for the master integrals. 

9. At four places information about the topologies are needed: as input for exp (cf. 
step 3) and crusher (cf. step 5), for the topology-specific FORM file (cf. step 4) 
and for identifying identical master integrals (cf. step 6). This input is generated 
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automatically from a file containing the definition of the three-loop topologies for 
the two-point on-shell integrals entering, e.g., the MS-on-shell relation of the quark 
mass [14-17]. This file is available from our previous calculation [17]. Note that 
this is the only input which contains non-trivial problem-specific information; the 
remaining input files (see steps 1 and 2) are either quite generic (e.g. identical for all 
QCD processes) or quite simple to adapt like the specification of the process under 
consideration for QGRAF. 

Let us stress that the automated setup outlined above requires little interaction from 
outside and thus minimizes possible errors. 

Many steps of the above list have been used extensively in previous calculations, 
however, the automatic calculation of the master integrals is new. Such an approach 
is essential in those cases where many master integrals occur. In our calculation we 
encounter 24 master integrals for the contribution involving a closed heavy fermion loop. 

Of course, a setup as described above requires several checks to be performed on the 
final result. On one hand they certainly include gauge parameter independence and the 
finiteness, on the other hand several checks on the numerical stability of our result are 
necessary. In our case the latter include the following: 

• FIESTA [13] , which is an efficient implementation of the sector decomposition method 
for the evaluation of master integrals, provides the possibility to introduce a lower 
cut-off in the numerical integration where numerical instabilities can occur. Actu- 
ally, for the evaluation of our master integrals we have to choose a non-zero cut-off. 
Its variation provides an estimate of the uncertainty. 

• FIESTA furthermore provides an uncertainty from the underlying Monte-Carlo inte- 
gration performed with Vegas [18] which we also take into account. 

• A further estimate of the uncertainty is provided by changing the basis used for the 
master integrals. The corresponding relations among the master integrals can be 
obtained in a straightforward way from the reduction tables generated by crusher. 

• Some of the master integrals are known analytically and can thus be used to replace 
the corresponding numerical expressions. 

3 Matching coefficient 

Starting from Eq. ([T]) it is possible to derive the equation^ 

Z 2 T V = c v Z v 1 , (3) 

where T v denotes the one-particle irreducible vertex diagrams with on-shell quarks with 
momenta q\ and q 2 and q 2 = (q± + q 2 ) 2 = 4mg. Some sample Feynman diagrams con- 
tributing at the one-, two- and three-loop level are shown in Fig.[IJ Z 2 is the wave function 

1 See, e.g., Ref. [7] for a derivation based on the threshold expansion [19,20]. 
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Figure 1: Feynman diagrams contributing to the matching coefficient. Solid lines de- 
note (heavy) quarks with mass itlq and curly lines denote gluons. In (d) and (h) mass 
counterterm diagrams are shown. 



renormalization constant in the on-shell scheme and Z v collects the infra-red divergences 
within the minimal subtraction scheme. The latter are cancelled against ultra-violet di- 
vergences of the effective theory rendering physical quantities finite. 

The two- and three-loop corrections to Z2 have been computed in Refs. [21] and [17,22], 
respectively, and Z v can be obtained from Ref. [23] (see also Refs. [7,24]). The latter reads 



Z„ = 1 



TT 



±fl + \c F cA 



TT 



h° F + h CA ) ^ 



25 



37 



C F + — C A — 

324 F 432 A e 



TT 



+ C F n h - > + 



(4) 



where C A = N c , C F = (iV c 2 - 1)/(2N C ) and T = 1/2 for a SU(N C ) group and the ellipses 
stand for non-fermionic and O (af ) terms. Note that the strong coupling is defined in the 
effective theory with n\ active quarks where n\ + is the total number of quark flavours. 
In our case we have n h = 1, however, we keep in the formulae for convenience. Since 
there are poles starting from order a 2 s1 higher order terms in e are necessary for the 
decoupling relation of a s . They have been computed in Ref. [25] and can be found in 
explicit form in Eq. (12) of Ref. [26]. For the evaluation of r„ to three- loop order also 
one- and two-loop expressions for the strong-coupling and quark-mass counterterms are 
needed which have been well-known for many years (see, e.g., Ref. [17]). 

The approach which has been chosen for the evaluation of the Feynman diagrams both 
at two-loop [27, 28] and at three-loop order [7] for the light-fermion contribution is based 



5 



on a partial fractioning of the integrands. As a consequence the occurring integrals can 
be mapped to diagrams containing less lines. Although they are in general simpler to 
evaluate they contain denominators raised to higher powers which partially compensates 
this advantage. On the other hand, the original vertex diagrams are more complicated 
to evaluate, but we can use additional recurrence relations derived from the fact that 
the momenta of the external quarks are the same. Actually, we computed the three- 
loop ^-contribution in both ways and observed no significant difference in the overall 
performance. Thus, in the automatic approach we refrain from performing the partial 
fractioning and evaluate directly the vertex diagrams. 

It is convenient to cast the perturbative expansion of the matching coefficient in the 
form 

= 1+ ^ 41 ) + f^y 4 2, + ('^wy c (3, +0(a 4 )i (5) 

7T \ 7T J \ IT I 

(3) 

where we further decompose c v according to the colour structures as 

4 3) = C F Tni (C F cffl + C A c FA L + Tn h c F hl + Tni c F ll) 
+ C F Tn h (C F cffh + C A c F ah + Tn h c F hh) 

+ non-fermionic and singlet terms . (6) 

The one- [29] and two-loop [27, 28, 30] terms have been known for more than ten years. 
More recently also the three-loop corrections proportional to ni became available [7]. The 
corresponding results reac{§ 

c« = ~2C F , 

c, 



P = f-^ + ^-^ln2-^ C (3)W 
V 72 144 6 4 y 



+ —C F Trn - - 
18 2 



c f l 



cffl = 46.7(1) +(-— + — 7r 2 --7r 2 ln2 + -C(3))L M + — tt 2 LJ , 



17 61 9 2 o, 1 , 1 

1 7r : 

12 36 3 

,181 67 9 5 9 13 M/ A /ll I „ , , 

' ' » = * M " ' ' + I ~64 ~ 432 + 9" k2 + T C(3) L " + T + T2* ) L J' 



557 26 o / 44 4 2 , , 

Cfhl = 1 ^ + 1 K L 

162 81 V 27 27 ' 

163 4 2 H r 2 

CFLL " "162 "27^ "27 L ^9^' (7) 



2 Note that in Ref. [7] the result has been expressed in terms of the six-flavour coupling whereas here 

(5) 

we use a s ■ This explains the difference in the logarithmic part of the coefficient cfhl- 
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this work 


Ref. [7] 


~{2) 
Cv 


-42.5138(2) 


-42.5140 


CFFL 


46.692(1) 


46.7(1) 


cfal 


39.623(1) 


39.6(1) 


cfhl 


-0.27029(4) 


-0.27025 


cfll 


-2.46833(3) 


-2.46834 


g(3)| 


n { (120.660(3) - 0.8228^) 


nj (121. - 0.8228^) 



Table 1: Two- loop and light-fermion contribution to c v . The results obtained with our 
numerical approach are compared to the ones of Ref. [7]. We have chosen the value 10~ 3 
for the FIESTA parameter If Cut. is set to one in the last row. 



with L M = ln^/m^) and (3q = (11/3 CU — 4/3Tnj)/4. In this paper we consider the 
contributions proportional to n^, i.e. new results for cffh, cfah and cfhh are presented. 
We separate the logarithmic contributions and write 



x 2 r 

Cffh — cffh — , 



1 

20 



I 121 11 2 , , 

Cfah = Cfah + - — vr L 



27 27 



cfhh — cfhh 



where the L M term in cffh arises from Eq. (TjO) and the one in cfah originates from the 
running of a s . 

Before considering the heavy-fermion contribution let us in a first step consider the two- 
loop and light-fermion contribution in order to get some confidence in our approach. In 
Tab. [T] the results of our approach are compared to the ones of Refs. [7]. The uncertainties 
given in the middle column are obtained by adding the numerical uncertainties of the 
individual master integrals in quadrature. We observe an impressive agreement for all 
coefficients which is in particular true for the analytically known coefficients c v , cfhl 
and cfll- Furthermore, in the case of £ffl and £fal a more precise result is obtained as 
compared to the approach of Ref. [7] where the Mellin-Barnes method has been used for 
the evaluation of the non-trivial master integrals. 

The new results for the contribution are shown in Tab. [2j In the column "numerical 
results" we show the numbers as obtained from the setup described in the previous section, 
i.e. there has been no manual interaction in the calculation of three-loop diagrams. 
On the other hand, if the master integrals which are known analyticallj{§ are used we 
obtain the numbers presented in the third column of Tab. [2J One observes only marginal 
improvements. For this reason we use for the following discussion the numerical results 
for the master integrals. 



3 The numerical agreement between the analytically known integrals and the results from FIESTA 
is about four to five digits in the coefficient of the highest e expansion term which enters the finite 
contribution of c v . 
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numerical 
result 



semi-analytical 
result 



CFFH"|log 
CFFH 



cfhh 



cfah 



cfhl 



-0.496(2) 
-0.841(3) 
-0.10(2) 
0.05126(1) 
-0.27029(4) 



-0.494(1) 
-0.840(2) 
-0.09(2) 
0.05124 
-0.27025 




0.93(4) -0.09010(l)n/ 



0.92(4) - 0.09008^ 



Table 2: Three-loop heavy fermion contribution to c v . For the results in the middle column 
we use only numerical results for the master integrals whereas in the right columns all 
available analytical information is employed. We have chosen the value 10~ 3 for the 
FIESTA parameter If Cut. is set to one in the last row. 



We have performed several checks on the correctness of our result and on the stability 
of the numerical calculations. In the Feynman rule for the gluon propagator we allow for 
a general gauge parameter and perform an expansion up to the linear term before the 
reduction to master integrals. The cancellation of the gauge parameter in the final result, 
after including counterterm contributions, serves as a welcome check for the correctness 
of our result. Furthermore, we have checked that both the spurious poles of order 1/e 4 
and 1/e 5 and the 1/e 3 and 1/e 2 poles in our final expression cancel with an accuracy of 
about 10 -4 . 

As a further test on the numerical stability of the evaluation of the master integrals 
we vary the parameter If Cut. In Fig. [2] the results for cffh and cfah are shown for 
several values between If Cut= 5 • 10~ 5 and If Cut= 0.1 where for the guidance of the eye 
the data points are connected by straight lines. One observes a broad plateau with only 
very minor variations. Larger deviations are obtained at the end points where either the 
value for If Cut becomes too big or numerical instabilities occur at the lower end of the 
integration region. Note that the CPU time for the evaluation of the master integrals 
with FIESTA varies from one to several days, depending on the setting for Vegas. 

As a further check on the numerical evaluation of the master integrals we have used 
a different momentum assignment in the input for FIESTA. As a consequence different 
expressions are generated in intermediate steps leading to different numerical integrations. 
The final results are in agreement with the ones in Tab. [2] within a one sigma level. 

A strong check on the numerical results for the master integrals is provided by a change 
of the master integral basis. We replace the complicated integrals by other integrals which 
we again evaluate with FIESTA. The relations between the old and new integrals can be 
extracted from the tables produced by crusher. By changing the basis two times (cf. 
Eqs. f[T3"j) and ffT4"|) ) we obtain the results given in Tab. [3J We find good agreement with 
the results of Tab. [2] for all coefficients. Let us mention that in the standard basis spurious 
poles of at most fourth order arise whereas both for basis 2 and 3 1/e 5 poles are present in 
our result. Let us furthermore stress that the integrals for the new bases are significantly 
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Figure 2: Coefficients cffh and cfah as a function of IfCut (choosing 
VarExpansionDegree=l). The band corresponds to the numerical uncertainty as pro- 
vided by FIESTA. 

more complicated due to the higher powers of the propagators. This explains the slightly 
worse numerical precision of the results in Tab. [3j 



Our final result for reads 

cffh 
cfah 



cfhh 



-0.841(6) , 
-0.10(4), 
427 158 2 16 

+ 7T 2 + — C 3) 

162 2835 9 w 



0.05124, 



(9) 



where the values for cffh and cfah are taken from Tab. [2j The uncertainties are con- 
servatively estimated by doubling the error from the numerical integration. In this way 
we account for effects connected to the FIESTA parameter If Cut, to different momenta 
assignments in the input of FIESTA, and the change of the master integral basis. The need 
for doubling the error can also be seen by comparing the two values for cfhh in Tab. [2j 
Inserting the numerical values for the colour factors we obtain for fi = itlq 



.(3) 



-0.823 nf + 120.66(1) m - 0.93(8) + non-fermionic and singlet terms . (10) 
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basis 2 



basis 3 



Eg. m 



Eg. (tHj) 



CFFHllog 
CFFH 



CFHH 



CFHL 



CFAH 



-0.50(1) 
-0.85(6) 
-0.15(9) 
0.0513(1) 
-0.27028(3) 



-0.496(8) 
-0.86(2) 
-0.13(4) 
0.0513(1) 
-0.2703(2) 




1.04(23) -0.09009(1)^ 



1.00(11) - 0.09009(5)71! 



Table 3: Three-loop heavy fermion contribution to c v . For the master integrals which are 
only available in numerical form the sets given in Egs. ({TBI and ffT4j) have been used. We 
have chosen the value 10 -3 for the FIESTA parameter If Cut. rih is set to one in the last 
row. 



It turns out the numerical coefficient of the heavy-fermion contribution is comparable 
with the nf part, however, significantly smaller than the coefficient of the linear ni term. 

4 Conclusions and outlook 

In this paper we discussed a completely automated approach for the calculation of three- 
loop vertex corrections contributing to the matching coefficients of the vector current, c v . 
In particular we consider the Feynman diagrams containing a closed heavy fermion loop 
which lead to significantly more complicated integrals than the light-fermion contributions 
considered more than two years ago [7] . We have shown that numerically stable results are 
obtained even for the case where all e coefficients of all 24 master integrals are evaluated 
numerically. Furthermore, we were able to improve the precision of the light-fermion 
contribution. The method developed in this work will be crucial for the three-loop non- 
fermionic and singlet contributions to c v which are still unknown. 

Our automated approach depends crucially on the fact that FIESTA is able to evaluate 
the master integrals to a sufficiently high accuracy. To make sure of this, we have checked 
that FIESTA reproduces all analytically known results within the given errors. In addition, 
we have performed the calculation of the matching coefficient in three different master 
integral bases and find consistent results. This makes us confident that our result is 
correct within the given error bar. 
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np2b2 npl2 np33 

Figure 3: Generic vertex diagrams where the solid and dashed lines correspond to massive 
and massless propagators, respectively. The external momenta are set to 4mg (out-going 
on the right-hand side) and toq (in-coming on the left-hand side), respectively. For V np \ 2 
and V np 33 closed fermion contributions are generated by contracting massless lines. 

Appendix: master integrals 

In this appendix we list the 24 master integrals which we encounter in the calculation for 
the n b contribution to c v . We refrain from providing explicit results since they can either 
be found in analytical form in the cited literature or can be obtained in numerical form 
using FIESTA [13]. 

There are twelve master integrals which are known in analytical form to a sufficiently 
high power in the e expansion. They are given by 

V np2b2 (0, 1, 0, 1, 1, 0, 0, 1, 1) , V np2b2 (0, 1, 0, 1, 1, 1, 0, 1, 1) , V np2b2 (0, 1, 1, 0, 1, 1, 0, 0, 0) , 
V np2b2 (0, 1, 1, 0, 1, 1, 1, 0, 1) , V np2b2 (0, 1, 1, 1, 1, 0, 0, 0, 0) , V np2b2 (0, 1, 1, 1, 1, 0, 0, 1, 1) , 
V np2b2 (0, 1, 1, 1, 1, 1, 0, 0, 0) , V npl2 (0, 1, 1, 1, 0, 0, 0, 0, 0) , v npl2 (o, i, i, o, 1, 0, 0, 0, 1) , 
V npl2 (0, 1, 1, 0, 1, 1, 0, 0, 0) , V npl2 (0, 1, 1, 1, 0, 1, 0, 0, 0) , V npl2 (l, 0, 1, 1, 0, 0, 0, 1, 0) , 

(11) 

where V np2b2 and V np i 2 are given in graphical form in Fig. [31 The numbers next to the 
lines mark the corresponding indices. The results for these integrals can be found in 
Refs. [22,33] (see also [7]). 

The poles of the integral V np2b2 (l, 0, 1, 0, 1, 1, 0, 0, 0) can be found in Ref. [34]; the finite 
part is computed with FIESTA. 

The following eleven integrals are only known numerically with the help of FIESTA 

V np2b2 (l, 0, 1, 0, 1, 1, 1, 0, 1) , V np2b2 (l, 0, 1, 0, 1, 2, 0, 0, 0) , V np2b2 (l, 0, 1, 1, 0, 1, 0, 1, 0) , 
V np2b2 (l, 0, 1, 1, 0, 1, 0, 2, 0) , V np2b2 (l, 0, 1, 1, 0, 1, 1, 1, 0) , V np2b2 (l, 0, 1, 1, 1, 1, 0, 0, 0) , 
V np2b2 (l : 0, 1, 1, 1, 1, 0, 1, 0) , V np2b2 (l, 0, 1, 1, 1, 1, 0, 2, 0) , k p33 (o, 1, 1, 1, 1, 1, 1, 0, 0) , 

K p33 (i, i, o, i, i, i, i, o, o) , k p33 (i, i, i, o, i, i, o, o, i) . (12) 

The change of the master integral basis discussed in Section [3] affects only the integrals 
in Eq. ( 1121) . V np2b2 (l, 0, 1, 0, 1, 1, 0, 0, 0) and the last five integrals in Eq. ( flTT) which are 



11 



replaced by either ( "basis 2" ) 

V npl2 (0, 2, 1, 0, 1, 0, 0, 0, 1) , V npl2 {0, 2, 1, 0, 1, 1, 0, 0, 0) , V npl2 (0, 2, 1, 1, 0, 1, 0, 0, 0) , 
V npl2 (l, 0, 2, 1, 0, 0, 0, 1, 0) , V npl2 (2, 0, 1, 1, 0, 0, 0, 1, 0) , V np2b2 (l, 0, 1, 0, 1, 1, 1, 0, 2) , 
V np2b2 (l, 0, 1, 0, 1, 2, 1, 0, 1) , V np2b2 (l, 0, 1, 0, 1, 3, 0, 0, 0) , V np2b2 (l, 0, 1, 1, 0, 1, 0, 3, 0) , 
V np2b2 (l, 0, 1, 1, 0, 1, 1, 2, 0) , V np2b2 {l, 0, 1, 1, 0, 1, 2, 1, 0) , V np2b2 (l, 0, 1, 1, 1, 1, 0, 3, 0) , 
V np2b2 (l, 0, 1, 1, 1, 2, 0, 0, 0) , V np2b2 (l, 0, 1, 1, 1, 2, 0, 1, 0) , V np2b2 (l, 0, 1, 1, 1, 2, 0, 2, 0) , 
V np33 (0 } 2, 1, 1, 1, 1, 1, 0, 0) , V np33 (2, 1, 0, 1, 1, 1, 1, 0, 0) , V np33 (2, 1, 1, 0, 1, 1, 0, 0, 1) , 

(13) 

or ("basis 3") 

V npl2 (0, 1, 1, 0, 1, 0, 0, 0, 2) , V npi3 (p, 1, 1, 0, 1, 2, 0, 0, 0) , V npl2 {0, 1, 1, 1, 0, 2, 0, 0, 0) , 
V npl2 (l, 0, 1, 1, 0, 0, 0, 2, 0) , V npl2 (l, 0, 2, 1, 0, 0, 0, 1, 0) , V np2b2 (l, 0, 1, 0, 1, 1, 2, 0, 1) , 
V np2b2 (l, 0, 1, 0, 2, 1, 0, 0, 0) , V np2b2 (l, 0, 1, 0, 2, 1, 1, 0, 1) , V np2b2 (l, 0, 1, 1, 0, 2, 0, 1, 0) , 
V np2b2 (l, 0, 1, 1, 0, 2, 0, 2, 0) , V np2b2 (l, 0, 1, 1, 0, 2, 1, 1, 0) , V np2b2 (l, 0, 1, 1, 1, 2, 0, 1, 0) , 
V np2b2 (l, 0, 1, 1, 2, 1, 0, 0, 0) , V np2b2 (l, 0, 1, 1, 2, 1, 0, 1, 0) , V np2b2 (l, 0, 1, 1, 2, 1, 0, 2, 0) , 
K P 3 3 (0, 1, 1, 1, 1, 1, 2, 0, 0) , V np33 (l, 1, 0, 1, 1, 1, 2, 0, 0) , V np33 (l, 1, 1, 0, 1, 1, 0, 0, 2) . 

(14) 

Note that both basis 2 and 3 contain one more master integral than our standard basis. 
This is because the latter in principle also contains the integral V np2b2 (l, 0, 1, 1, 1, 2, 0, 1, 0), 
however, for our application the corresponding coefficient is zero. 

References 



[i] 

[2] 

[3] 
[4] 
[5] 

[6] 

[7] 



W. E. Caswell and G. P. Lepage, Phys. Lett. B 167 (1986) 437. 

G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51 (1995) 1125 [Erratum- 
ibid. D 55 (1997) 5853] |arXiv:hep-ph/9407339| . 



N. Brambilla et al, arXiv:hep-ph/0412158 



M. Martinez and R. Miquel, Eur. Phys. J. C 27 (2003) 49 [arXiv:hep-ph/0207315] . 

M. Beneke, Y. Kiyo and K. Schuller, PoS RADCOR2007 (2007) 051 
jarXiv:0801.3464l [hep-ph]]. 

A. V. Smirnov, V. A. Smirnov and M. Steinhauser, Phys. Lett. B 668 (2008) 293 
jarXiv:0809. 19271 [hep-ph]]. 



P. Marquard, J. H. Piclum, D. Seidel and M. Steinhauser, Nucl. Phys. B 758 (2006) 



144 arXiv:hep-ph/0607168 



12 



P. Nogueira, J. Comput. Phys. 105 (1993) 279. 



[9] J. A. M. Vermaseren, arXiv:math-ph/0010025 



[10] R. Harland er, T. Seidensticker and M. Steinhauser, Phys. Lett. B 426 (1998) 125 
[hep-ph/9712228| . 



[11] T. Seidensticker, hep-ph/9905298. 



[12] P. Marquard and D. Seidel, unpublished. 

[13] A. V. Smirnov and M. N. Tentyukov, larXiv:0807.41~29l [hep-ph]. 

[14] K. G. Chetyrkin an d M. Steinhauser, Nucl. Phys. B 573 (2000) 617 
[arXiv:hep-ph/9911434| . 

[15] K. G. Chetyrkin an d M. Steinhauser, Phys. Rev. Lett. 83 (1999) 4001 
[arXiv:hep-ph/9907509l . 

[16] K. Melnikov and T. v. Ritbergen, Phys. Lett. B 482 (2000) 99 
[arXiv:hep-ph/99123"91~l . 

[17] P. Marq uard, L. Mihaila, J. H. Piclum and M. Steinhauser, Nucl. Phys. B 773 (2007) 
1 [arXiv:hep-ph/0702185] . 

[18] G. P. Lepage, J. Comput. Phys. 27 (1978) 192. 



[19] M. Beneke and V. A. Smirnov, Nucl. Phys. B 522 (1998) 321 |arXiv:hep-ph/9711391 



[20] V. A. Smirnov, "Applied asymptotic expansions in momenta and masses," Springer 

(2002) . 

[21] D. J. Broadhurst, N. Gray and K. Schilcher, Z. Phys. C 52 (1991) 111. 

[22] K. Melnikov and T. van Ritbergen, Nucl. Phys. B 591 (2000) 515 
[arXiv:hep-ph/0005131| . 

[23] B. A. Kniehl, A. A. Penin, M. Steinhauser and V. A. Smirnov, Phys. Rev. Lett. 90 

(2003) 212001 [Erratum-ibid. 91 (2003) 139903(E)] [arXiv:hep-ph/0210161] . 

[24] M. Beneke, Y. Kiyo and A. A. Penin, Phys. Lett. B 653 (2007) 53 [ arXiv: 0706 .27331 
[hep-ph]]. 

[25] K. G. Chetyrkin, B. A . Kniehl and M. Steinhauser, Nucl. Phys. B 510 (1998) 61 
[arXiv:hep-ph/9708255| . 

[26] A. G. Grozin, P. Marquard, J. H. Piclum and M. Steinhauser, Nucl. Phys. B 789 
(2008) 277 |arXiv:0707. 13881 [hep-ph]]. 



13 



A. Czarnecki and K. Melnikov, Phys. Rev. Lett. 80 (1998) 2531 
|arXiv:hep-ph/9712222| . 

M. Beneke, A. Signer and V. A. Smirnov, Phys. Rev. Lett. 80 (1998) 2535 
|arXiv:hep-ph/9712302| . 

G. Kallen and A. Sarby, K. Dan. Vidensk. Selsk. Mat.-Fis. Medd. 29, N17 (1955) 1. 

B. A. Kniehl, A. Onishchenko, J. H. Piclum and M. Steinhauser, Phys. Lett. B 638 



(2006) 209 [arXiv:hep-ph70 604072|. 

J. A. M. Vermaseren, Comput. Phys. Commun. 83 (1994) 45. 

D. Binosi and L. Theussl, Comput. Phys. Commun. 161 (2004) 76 
|arXiv:hep-ph/0309015| . 

S. Laporta and E. Remiddi, Phys. Lett. B 379 (1996) 283 |arXiv:hep-ph/9602417| . 

S. Groote, J. G. Korn er and A. A. Pivovarov, Eur. Phys. J. C 36 (2004) 471 
|arXiv:hep-ph/0403122| . 



14 



